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Abstract 

The  paper  considers  the  development  and  simulation 
testing  of  the  control  algorithms  for  an  autonomous 
high-glide  aerial  delivery  system,  which  consists  of 
650 sq.ft  rectangular  double- skin  parafoil  and  500 Ih 
payload.  The  paper  applies  optimal  control  analysis  to 
the  real-time  trajectory  generation.  Resulting  guidance 
and  control  system  includes  tracking  this  reference  tra¬ 
jectory  using  a  nonlinear  tracking  controller.  The  paper 
presents  the  description  of  the  algorithm  along  with 
simulation  results  and  ends  with  guidelines  for  the  fur¬ 
ther  implementation  of  the  developed  software  aboard 
the  real  system. 

I.  Introduction 

Maneuverable  ram-air  parafoils1'3  are  widely  used  to¬ 
day.  The  list  of  users  includes  sky  divers,  smoke  jump¬ 
ers  and  special  forces.  Furthermore,  their  extended 
range  (compared  to  that  of  round  canopy  parachutes) 
makes  them  very  practical  for  payload  delivery.  Recent 
introduction  of  the  Global  Positioning  System  (GPS) 
made  the  development  of  fully  autonomous  ram-air 
parafoils  possible.  Moreover,  ram-air  parafoils  are 
nowadays  considered  a  vital  element  of  the  unmanned 
air  vehicles  (UAVs)  capability  and  of  International 
Space  Station  (ISS)  crew  rescue  vehicle.4'12 

Autonomous  parafoil  capability  implies  delivering 
the  system  to  a  desired  landing  point  from  an  arbitrary 
release  point  using  onboard  computer,  sensors  and  ac¬ 
tuators.  This  requires  development  of  a  guidance,  navi¬ 
gation  and  control  (GNC)  system.  The  navigation  sub¬ 
system  manages  data  acquisition,  processes  sensor  data 
and  provides  guidance  and  control  subsystems  with 
information  about  parafoil  states.  Using  this  informa¬ 
tion  along  with  other  available  system  data  (such  as 
local  wind  profiles),  the  guidance  subsystem  plans  the 
mission  and  generates  feasible  (physically  realizable 
and  mission  compatible)  trajectory  that  takes  the  para- 


This  paper  is  declared  a  work  of  the  U.S.  Government  and  is 
not  subject  to  copyright  protection  in  the  United  States. 

Associate  Professor;  Senior  Member  AIAA. 

^Research  Associate  Professor,  Associate  Fellow  AIAA. 


foil  from  the  initial  position  to  the  desired  landing 
point.  Finally,  it  is  the  responsibility  of  the  control  sys¬ 
tem  to  track  this  trajectory  using  the  information  pro¬ 
vided  by  the  navigation  subsystem  and  onboard  actua¬ 
tors. 

In  the  past  decade,  several  GNC  concepts  for  glid¬ 
ing  parachute  applications  have  been  developed  and 
published.13'22  Most  of  them  were  tested  in  a  simulation 
environment,  some  in  flight  test. 

Some  of  the  published  papers  include  analytical 
development  of  feasible  trajectories.  In  Ref. 23,  for  ex¬ 
ample,  minimum  glide  angle  and  descent  rate,  mini¬ 
mum  turn  radius  and  maximum  gliding  turn  rate,  reach¬ 
able  ground  area,  minimum  altitude  loss  problems  were 
considered  and  solved  analytically  for  a  simplified 
three-degree-of-freedom  (3DoF)  model  with  a  parabolic 
lift-drag  polar.  The  problems  of  the  maximum  time  to 
stall  the  glider  at  the  predefined  point  in  a  horizontal 
plane  and  a  maximum  range  for  the  appropriate  simpli¬ 
fied  model  with  constraints  were  also  considered.  Op¬ 
timization  included  no  wind  scenario. 

Assuming  the  wind  velocity  to  be  constant  Ref.24 
synthesized  zero-lag  control  of  the  gliding  parachute 
system  with  the  constant  rate  of  descent.  This  interest¬ 
ing  work  uses  a  compound  performance  index  blending 
time,  kinetic  energy,  and  controls  expenditure.  The  op¬ 
timal  problem  was  solved  by  applying  Pontrjagin's 
maximum  principle.25  The  resulting  two-point  bound¬ 
ary-value  problem  was  solved  using  Krylov- 
Chernous'ko  method.26 

However,  majority  of  the  published  results  on 
GNC  system  development  for  parafoils  sue  heuristic 
rules  and  classical  control  algorithms.  It  is  appropriate 
to  mention  some  of  them  here  emphasizing  a  scale  of 
the  system,  sensor  suite  specs  and  control  approach. 

Probably  one  of  the  first  GNC  algorithms  for  small 
ram- air  parachutes  is  presented  in  Ref.  13.  Computer 
programs  were  written  using  6DoF  models  of  small 
parafoils,  such  as  Ram-air,  Para-commander,  and  Stan¬ 
dard  Round.  The  aerodynamic  models  were  simplified, 
but  they  included  the  non-linearities  and  severe  con¬ 
straints  common  to  all  parachutes  as  compared  to  air¬ 
craft.  The  algorithms  used  GPS  position  and  velocity 
data  provided  at  1  Hz  to  control  both  toggle  lines.  Hom¬ 
ing  and  circling  algorithms  were  simulated  using  direc- 
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tion  sines  and  cosines.  Wind  was  included  in  all  algo¬ 
rithms.  Homing  algorithms  were  included  that  in¬ 
creased  the  allowable  offset,  minimized  the  miss  dis¬ 
tance  and  reduced  the  impact  velocity.  These  were 
based  on  two-axis  proportional  navigation  (PN),  and 
also  upon  upward-biased  rate- integrating  navigation  in 
glide  angle.  The  latest  produced  final  flare  and  nearly 
vertical  landing  even  from  a  long  initial  offset.  For  the 
case  of  high  initial  altitude  and  small  offset  circling 
algorithms  were  also  introduced.  Finally,  combined 
downwind,  upwind,  and  guidance  algorithms  were 
simulated  to  emulate  manned  parachute  procedures. 

Another  approach  was  employed  in  Ref.16  where 
Pseudo-Jumper  (PJ)  guidance  and  control  algorithms 
based  on  the  professional  jump  procedures  were  intro¬ 
duced.  The  main  phases  of  the  mission  after  release 
were  defined  as  follows: 

•/"Cruise"  phase  consisting  of  aiming  for  the  holding 
area  way-points  (W  Ps)  by  one  or  more  types  of  hom¬ 
ing  algorithms; 

/"Hold"  phase:  As  the  parachute  gets  close  enough  to 
the  homing  WPs,  it  is  caused  to  perform  S-turns  just 
downwind  of  the  approach  zone.  The  pattern  is  ad¬ 
justed  to  suit  the  wind  speed  so  that  it  stays  near  the 
landing  zone  (LZ)  and  has  just  enough  upwind  dis¬ 
tance  left  for  the  approach  and  flare; 

/"Center"  phase:  at  some  point,  governed  mainly  by 
altitude  (time-to-go),  the  PJ  must  aim  for  the  center 
of  the  drop  zone  so  that  it  is  heading  towards  the 
touchdown  target.  Enough  time  must  be  allowed  to 
complete  any  ongoing  S-turn; 

/"Approach"  phase:  at  some  point,  the  PJ  has  to  stop 
trying  to  center  and  get  lined  up  with  the  surface 
wind  direction  and  to  minimize  toggle  actions  past 
that  point; 

/"Flare"  phase  is  triggered  with  rapid  full  braking  to 
arrest  the  vertical  and/or  horizontal  inertial  velocities 
in  an  optimal  control  fashion. 

Perfect  knowledge  of  parafoil's  relative  location, 
i.e.  availability  of  a  GPS  sensor  and  of  sonic  altimeter 
as  well  as  of  winds  below  300ffwas  assumed. 

Recent  work  by  the  Iowa  State  University 
(http://cosmos.ssol.iastate.edu/RGS)27,28  and  Atair 
A  erospace  (http://www.extremefly.com/aerospace) 
includes  development  of  a  parafoil-based  Recovery 
Guidance  System  to  return  payloads  from  a  high- 
altitude  balloon  back  to  the  ground  safely  and  predicta¬ 
bly.  The  parafoil  is  supposed  to  accompany  the  balloon 
and  to  be  cut  free  of  the  balloon  at  a  predetermined 
time.  It  then  guides  the  payload  to  a  designated  LZ.  As 
reported,  fully  automated  guidance  hardware  and  soft¬ 
ware  used  DGPS  (differential  GPS)  data,  a  digital  com¬ 
pass,  a  linear  quadratic  regulator  and  PN  guidance  to 
control  and  navigate  the  parafoil.  Only  simulation  re¬ 
sults  have  been  reported  so  far. 


Another  GNC  system  incorporating  appropriate 
data  acquisition  system  was  developed  to  support  high- 
altitude-high-opening  insertion  into  hostile  environ¬ 
ments  and  was  supposed  to  be  used  with  a  500 lb  para¬ 
foil.  The  GPS  guidance  mode  was  designed  to  steer  a 
glider  to  the  LZ  using  only  GPS  data  or  GPS  plus  up  to 
eight  additional  sensors.  The  GPS  data  was  provided  at 
1Hz  and  inertial  sensors  at  200 Hz.  Control  steering  was 
accomplished  by  actuating  one  of  the  two  steering  con¬ 
trols,  which  allowed  four  flight  modes:  straight  flight, 
steer  left  or  right,  hard  turn  left  or  right,  and  flare.  No 
test  results  however  are  available  so  far. 

A  nother  project  was  reported  by  the  German  A  ero¬ 
space  Center.  It  aimed  to  identify  behavior  of  a  small 
(200 lbs)  and  extensively  instrumented  parafoil-based 
system  and  to  investigate  GNC-concepts  for  the 
autonomous  landing.29,30  Three  and  four  degree  of  free¬ 
dom  models  were  developed  and  parameters  of  these 
models  were  identified  during  this  study. 

The  world's  largest  parafoil  was  designed  for  the 
U.S.  Army  by  Pioneer  Aerospace  Corp.  This  was  a  part 
of  a  contract  undertaken  by  Pioneer  and  the  U  ,S.  A  rmy 
Soldier  Systems  Command  and  was  tested  in  2001- 
2002  at  the  U  .S.  Y  uma  Proving  Grounds,  AZ  (Y  PG)  as 
an  aid  for  ISS  "lifeboat".  The  resulting  parafoil  was 
ultimately  intended  to  evacuate  a  seven-person  crew 
from  the  station  in  case  of  an  emergency  (Fig. 1). 11-12 

The  1,600/Jb  (50  times  the  weight  of  a  typical  per¬ 
sonal  parachute)  7,500sg. ft  parafoil  wing  has  the  wing¬ 
span  wider  than  a  football  field  and  a  canopy  area  of 
more  than  one  and  a  half  times  that  of  the  wings  of  a 
747  jumbo  jet  and  is  designed  to  carry  weights  of  up  to 
42,000/bs,  the  heaviest  cargo  ever  transported  by  a 
parafoil.  The  system  is  able  to  move  anything  from 
supplies  to  armored  vehicles  into  strategic  areas.  Other 
applications  could  include  mentioned  recovery  of  high- 
value  assets  for  the  space  program,  and  humanitarian 
aid,  such  as  the  delivery  of  food,  supplies  and  medica¬ 
tions  to  war-torn  areas  like  Bosnia,  Sarajevo,  Afghani¬ 
stan,  and  Iraq. 

The  Pioneer  Aerospace  Corporation  is  an  industry 
leader  in  parachute  design.  They  have  provided  recov¬ 
ery  parachutes  for  NASA's  space  shuttle,  the  Air 
Force's  B-1A  Bomber  and  parachutes  for  other  NASA 
and  Department  of  Defense  programs.  Since  testing  on 
the  Guided  Parafoil  Aerial  Delivery  System  (GPADS) 
program  began  in  January  1994,  Pioneer  holds  records 
for  the  largest  parafoil,  highest  altitude  (22,000ft),  and 
heaviest  weight  (15,500/bs)  airdropped. 

The  GNC  system  for  GPADS7,8  was  developed  by 
the  European  Space  Agency.22,29-33  The  Parafoil  Tech¬ 
nology  Demonstrator  program  was  set  up  and  con¬ 
ducted  under  ESA  to  demonstrate  feasibility  of  guiding 
a  large-scale  parafoil  autonomously  to  a  predefined 
point  and  to  perform  a  flared  landing  within  a  specific 
range.  Another  objective  was  to  enable  further  investi- 
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supposed  to  be  completely  automated,  although  the 
crew  will  have  the  capability  to  switch  to  backup  sys¬ 
tems  and  steer  the  parafoil  manually,  if  necessary. 


gation  into  the  flight  dynamics  of  the  parafoil  systems. 
The  first  automatically-guided,  parafoil  based  descent 
and  landing  of  a  test  vehicle  was  carried  out  in  support 
of  ESA's  development  in  1997.  The  spacecraft  was 
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Figure  1  The  X  -38  vehicle  under  the  world's  largest  parafoil  (courtesy  of  N  A SA ) 


The  prototype  reference  trajectory  (RT)  consisted 
of  the  following  nominal  phases 
/Acquisition  -  trajectory  initialization  and  initial  turn; 
/Homing  -  travel  to  intercept  Energy  Management 
Circle  (EMC); 

/EMC  Entry,  EMC  and  EMC  Exit  -  the  energy  man¬ 
agement  circling  maneuver; 

/  Final  -  final  approach  to  target  on  desired  heading; 

The  GNC  system  employs  laser  or  radar  altimeter 
to  trigger  flare  and  is  required  to  steer  the  parafoil  to 
within  100m  of  a  preselected  desired  impact  point  (DIP) 
(so  far  the  accuracy  of  200m  was  achieved  on  a  smaller 
scaled  prototypes). 

To  support  the  development  of  GPADS  software 
and  to  evaluate  its  expected  performance  on  ram-air 
parafoils  with  payload  capacities  ranging  from  less  than 
200/bs  to  40,000/bs,  Draper  Laboratory  also  developed 
a  GNC  package.17'19  This  software  was  demonstrated 
using  NASA  Dryden  Flight  Research  Center  airdrop 
test  article  known  as  the  "Wedge  3"  which  utilizes  an 
88 sq.ft  ram-air  parafoil  and  carries  a  payload  of  175/bs. 
Same  GNC  algorithms  were  also  tested  in  simulation 
using  the  model  of  a  larger  3,600 sq.ft  guided  parafoil 
air  delivery  system  with  a  13,088/b  payload.  Upon  full 
deployment  of  the  parafoil  canopy  the  GNC  algorithm 
is  activated.  It  determines  the  desired  heading  angle  that 


will  lead  to  a  touchdown  at  LZ.  Seven  guidance  modes 
are  available: 

/  Automatic  trim; 

/  M  aneuver  to  W  P; 

/  Transition  to  hold; 

/  Holding  pattern  (to  maintain  turn  radius); 

/  M  aneuver  to  target; 

/  Landing  approach; 

/  Landing  flare. 

The  control  system  may  switch  between  any  of  these 
modes  sequentially  or  bypass  some  of  them. 

The  guidance  concept  presented  in  Ref.22  follows 
generally  the  same  scheme.  The  energy  management 
concept  is  adopted  motivated  by  the  way  novice  sky- 
divers  learn  the  landing  approach.16  The  trajectory  is 
planned  using  W  Ps.  The  guidance  works  as  follows: 
/Homing:  from  the  current  position  the  energy  man¬ 
agement  center  (EMC)  is  approached  if  possible.  If 
the  remaining  distance  is  not  long  enough  to  reach 
the  desired  impact  point,  the  WP  is  moved  in  the  di¬ 
rection  to  the  final  turn  point  (FTP); 

•/Energy  Management:  after  reaching  EMC,  surplus 
altitude  is  reduced  by  flying  S-patters  to  the  energy 
management  turn  points  (EM TP).  If  the  remaining 
distance  becomes  too  small  to  reach  DIP,  the  EM  TP 
is  shifted  along  the  energy  management  axis  towards 
the  EMC; 
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^Landing:  the  landing  phase  is  divided  into  4  sub¬ 
phases:  (1)  transition  into  landing  corridor,  (2)  ap¬ 
proach  of  final  turn  point  (FTP),  (3)  turn  into  wind 
and  (4)  flare. 

The  Backup  mode  was  also  developed  as  follows: 
if  altitude  and  subsequently  remaining  distance  are  too 
small  to  reach  the  landing  point  during  homing  phase, 
the  vehicle  keeps  heading  the  FTP  until  a  specified  de¬ 
cision  altitude  is  reached.  Then  the  vehicle  turns  against 
wind  independently  of  its  position  and  prepares  for 
landing. 

The  guidance  system  output  consists  of  a  heading 
command  computed  using  current  position  and  next 
WP.  For  wind  accommodation,  the  trajectory  is  planned 
in  the  wind  fixed  coordinate  system.21  Assuming  con¬ 
stant  horizontal  wind,  the  actual  position  is  shifted  to  a 
virtual  position  by  the  predicted  wind  drift,  the  vehicle 
will  be  affected  during  the  remaining  flight  until  touch 
down.  Within  this  wind  fixed  system  the  guidance  can 
be  done  as  if  there  is  no  wind.  Depending  on  the  sign  of 
the  heading  error,  either  the  left  or  the  right  actuator  is 
activated,  leading  to  a  left  or  a  right  turn.  A  saturation 
limited  PN  controller  computes  the  desired  actuator 
position.  For  instance,  with  an  amplification  of  4,  a 
commanded  turn  of  90°  first  activates  one  actuator  to 
full  stroke  until  the  heading  difference  goes  below  15° 
and  the  proportional  part  of  the  controller  takes  over. 

According  to  Ref. 22  after  selective  availability  has 
been  switched  off  (in  May  2000),  plain  GPS  data  with¬ 
out  differential  corrections  is  accurate  enough  for  tra¬ 
jectory  planning.  The  GPS  altitude,  with  errors  of  about 
±20 m,  is  combined  with  the  pressure  altitude.  Pleading 
information  is  supposed  to  be  computed  by  applying  a 
complementary  filter  to  blend  azimuth  and  yaw  rate  as 
described  in  Refs.14-15.  Simulation  results  also  show 
that  possible  wind  mispredictions  have  more  influence 
on  the  landing  accuracy,  than  incorrect  estimates  of 
model  parameters  and  measurement  errors.  However,  as 
long  as  the  errors  stay  in  the  normal  range  (position 
±10m,  altitude  ±5 m,  heading  ±10°)  the  GNC  algorithm 
is  robust  enough  to  land  the  vehicle  within  a  range  of 
±50m  around  the  predefined  landing  point,  provided 
that  the  wind  is  perfectly  known. 

As  seen  above  in  most  GNC  concepts,  the  guidance 
itself  is  divided  into  three  phases:  homing,  energy  man¬ 
agement  and  landing,  where  each  of  them  can  consist  of 
further  subphases.  Homing,  i.e.  flying  towards  the  LZ, 
and  landing  are  similar  in  almost  all  guidance  concepts. 
Energy  management  is  done  by  flying  a  certain  pattern 
close  to  the  landing  point,  until  some  criteria  is  satis¬ 
fied,  switching  to  the  landing  phase.  Most  differences 
among  the  realized  guidance  concepts  are  concerned 
with  the  shape  of  the  pattern  and  the  criteria  for  the 
landing  decision. 

Present  paper  addresses  the  problem  of  GNC  de¬ 
velopment  for  a  parafoil  system  as  follows.  First,  a  fea¬ 


sible  trajectory  is  generated  in  real-time  that  takes  the 
parafoil  from  initial  release  point  to  touch  down  at  a 
desired  impact  point.  It  is  assumed  that  only  the  direc¬ 
tion  of  the  wind  at  LZ  is  known.  The  trajectory  consists 
of  an  initial  glide,  spiral  descent  and  of  final  glide  and 
flare.  The  final  glide  and  flare  are  directed  into  the  wind 
at  LZ.  This  structure  is  motivated  by  optimal  control 
analysis  carried  out  on  a  simplified  model.  The  result¬ 
ing  trajectory  is  tracked  using  a  nonlinear  algorithm 
with  guaranteed  local  stability  and  performance  proper¬ 
ties  (see  Ref.34).  The  control  algorithm  converts  trajec¬ 
tory-tracking  errors  directly  into  control  actuator  com¬ 
mands.  Therefore,  only  GPS  position  and  velocity  are 
needed  to  implement  it. 

As  expected  the  basic  trajectory  structure  is  similar 
to  the  ones  reported  in  the  literature.  However,  as 
shown  in  Section  III,  a  single  smooth  inertial  trajectory 
is  generated  using  simple  optimization  and  is  tracked 
throughout  the  drop  by  the  same  control  algorithm.  This 
eliminates  the  need  for  multiple  modes,  extensive 
switching  logic  and  wind  information  throughout  the 
drop.  The  latter  is  made  possible  by  the  nonlinear  con¬ 
trol  algorithm  that  tracks  the  inertial  trajectory  directly 
and  treats  wind  as  a  disturbance. 

This  specific  delivery  system  considered  in  this  pa¬ 
per  is  called  Pegasus.  It  consists  of  650 sq.ft  span  rec¬ 
tangular  double-skin  parafoil  and  500 lb  payload  (Fig. 2) 
and  was  originally  manufactured  by  the  FXC  Corp.  It 
can  be  controlled  by  symmetric  and  differential  flap 
deflections  that  occupy  outer  four  (out  of  eight  on  each 
side)  cells  of  the  parafoil. 

This  paper  is  organized  as  follows.  Section  1 1  con¬ 
tains  a  brief  description  of  the  system  and  of  the  com¬ 
plete  6DoF  model  used  to  support  the  development  of 
the  GNC  algorithm.  Section  III  introduces  optimal  con¬ 
trol  strategy  for  a  simplified  parafoil  model  using 
Pontrjagin's  principle  of  optimality.  Section  IV  dis¬ 
cusses  real-time  trajectory  generation  for  the  control 
algorithm.  Section  V  discusses  the  development  of  the 
tracking  control  algorithm.  Section  VI  introduces  simu¬ 
lation  results  of  the  complete  guidance  and  control  sys¬ 
tem.  Finally,  Section  VII  contains  the  main  conclusions. 

II.  System  Architecture  and  Modeling 
To  support  the  development  of  the  GNC  algorithm,  a 
complete  6DoF  model  of  Pegasus  was  developed  and 
tuned  using  flight  test  data  provided  by  theYPG.  Fig. 3 
outlines  Pegasus  dimensions. 

The  description  of  the  6DoF  model  used  is  given  in 
Ref. 35.  This  model  matches  flight  test  data  and  is  char¬ 
acterized  by  the  following  integral  parameters:  the  av¬ 
erage  descent  rate  is  3. 7-3. 9 m/s,  glide  ratio  is  about  3.0, 
the  turn  rate  of  ~6°/s  corresponds  to  the  full  deflection 
of  one  flap.  An  interesting  feature  from  the  control 
standpoint  is  that  the  system  exhibits  almost  no  flare 
capability  (flaps  deflection  results  in  almost  no  change 
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in  the  descent  rate).  Therefore,  it  is  very  important  to  land  it  into  the  wind. 


Figure  2  Double-skin  650sg. ft  parafoil  with  500/b  payload 
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Figure  3  The  5001b  parafoil  geometry 
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III.  Optimal  Control  Synthesis 

Consider  the  following  kinematic  model  of  a  parafoil  in 
the  horizontal  plane.  Suppose  we  have  a  constant  glide 
ratio  and  by  pulling  risers  we  can  control  its  yaw  rate. 
M  athematically,  this  is  expressed  by  the  following  sim¬ 
plified  equations: 

X  =  V  cos^,  y  =  V  sin ^ ,  \j/  =  v,  (1) 


(2) 


where  ve  [-£;£]  is  the  only  control. 

The  Hamiltonian  for  the  system  (1)  for  a  time- 
optimal  control  can  now  be  written  as: 

H  =P,^(p„p,)(“S£ 

where  equations  for  adjoint  variables  px,  py  and  p 
are  given  by 

Px  =  0,  Py=  0, 

'  s\ni// 


P¥=v(pxlpy) 


Aon...  a  (3) 


-COSy/ 


The  optimal  control  for  the  time-minimum  problem 
now  is  given  by 

v  =  zsign(p¥).  (4) 


By  differentiating  the  expression  for  p¥  (3)  and 

combining  it  with  Hamiltonian  (2)  for  both  cases  when 
p¥  >  0  and  p¥  <  0  we  can  get  a  set  of  equations  for 

Pr: 

P¥+Z2p¥TS  =  0.  (5) 


This  differential  equation  gives  two  sinusoids 
(shifted  with  respect  to  abscise  axis  by  +H_1)  as  solu¬ 
tions  for  the  general  (non-singular)  case 

pv  =C1sin(St  +  C2)±S“1,  (6) 

where  Cj  and  C2  are  constants  defined  by  the  concrete 
boundary  conditions.  If  C^ET1  the  parafoil  model 
moves  along  a  descending  spiral.  It  takes  2/zS-1  sec¬ 
onds  to  make  a  full  turn  with  a  radius  of  \/£-1.  If 
Cj  =S_1  there  exists  a  possibility  of  singular  control. 
This  is  caused  by  the  fact  that  there  exists  a  point  in 
time  where  both  pv  and  p¥  are  zero  as  can  be  seen 

from  (6). 

Consider  singular  control  for  this  model.  By  defini¬ 
tion  it  means  that  pv=  0 .  For  the  time-optimal  prob¬ 
lem  from  the  Hamiltonian  (2)  and  third  equation  in  (3) 
(of  course  keeping  in  mind  the  first  two)  it  follows  that 
for  a  singular  control  case 

px  =1/ _1  COSy/  ,  Py  =i/'1siny/  ,  y/  =  const .  (7) 

Expressions  (7)  imply  that  singular  control  corre¬ 
sponds  to  motion  with  a  constant  heading  (v  =  0).  It 
may  not  however  be  realized.  Instead,  the  parafoil 
model  may  switch  from  right-handed  spiral  to  a  left- 
handed  one  or  vice  versa.  Planar  projections  of  possible 
trajectories  are  shown  on  Fig.4.  Concrete  boundary 
conditions  define  which  case  will  be  realized. 


To  summarize,  Figs.5,6  show  examples  of  three- 
dimensional  (3D)  time-optimal  trajectories  computed  in 
the  inertial  frame  for  the  gliding  system  (1).  As  sug¬ 
gested  by  the  Pontrjagin's  principle  of  optimality  they 
would  consist  of  helices  (spiral  descent)  and  straight 
descent  segments.  The  basic  difference  between  two 
trajectories  shown  on  Fig. 5  is  whether  parafoil  expends 


its  potential  energy  at  the  beginning  of  descent  or  at  the 
end.  It  may  depend  on  the  tactical  conditions  in  the  LZ, 
terrain,  or  some  other  factors.  In  a  general  case  parafoil 
could  spend  its  potential  energy  in  the  vicinity  of  some 
other  WP  differing  from  the  start  and  final  portions  of 
the  trajectory  as  itshown  in  Fig. 6.  In  any  case,  since  the 
control  system  cannot  meet  any  time  constraints  due  to 
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unavailability  of  thrust  on  a  parafoil  system,  the  com-  inertial  frame  and  are  time  independent, 
mon  feature  of  these  RTs  is  that  they  are  defined  in  the 


Figure  5  3D  representation  of  possible  inertial  RTs 


A  nother  critical  issue  is  that  in  order  to  meet  a  soft- 
landing  requirement  the  parafoil  at  the  LZ  must  align 
itself  into  the  wind. 

IV.  Real-Time  Trajectory  Generation 

The  optimal  control  analysis  of  the  previous  section 
motivated  the  basic  RT  structure  as  is  shown  in  Fig.  7. 
This  RT  consists  of  three  segments:  initial  straight-line 
glide  (segment  1),  spiral  descent  (segment  2)  and  final 
glide  and  flare  (segment  3).  Since  it  is  required  that  the 
parafoil  is  aligned  into  the  wind  at  LZ,  the  final  glide 
segment  ends  at  the  desired  impact  point  (DIP)  and  is 
directed  into  the  wind.  Furthermore,  the  final  glide 
starts  an  offset  distance  doffset  away  and  a  certain  height 

above  DIP.  The  height  is  determined  by  the  flight  path 
angle  of  the  segment  3.  Similarly,  the  first  segment 
starts  at  the  parafoil  release  point  (RP)  and  ends  at  a 
point  defined  by  the  flight  path  angle  of  the  first  seg¬ 
ment  at  the  spiral  descent  segment.  The  first  and  last 
segments  are  fused  together  by  the  spiral  descent  seg¬ 


ment.  The  radius  of  the  spiral  is  adjusted  to  provide 
smooth  transition  between  each  segment.  The  rest  of 
the  section  derives  the  mathematical  representation  of 
the  complete  RT. 

Let  yRTi  denote  the  desired  flight  path  angle  for  the 

first  straight-line  segment  of  the  trajectory,  similarly  let 
yRh  denote  the  flight  path  angle  for  the  spiral -descent 

segment  and,  finally,  let  yRT  represent  the  flight  path 

angle  for  the  final  straight-line  segment.  Furthermore, 
let  y/f  represent  the  wind  direction  at  the  LZ, 

Po  =(xo-y0-zo)r  is  the  RP,  and  p f  =(xf,yf,zf)T  de¬ 
notes  the  DIP.  Define  vector 
i  =  (cos(y/>  +7r),s\r\(y/f  +^),0)r  to  be  the  unit  vector 
that  points  into  the  wind  at  the  LZ  and  vector 
i*  =  (-sin(y/>  +  n), cos[y/f  +^),0)T  to  be  the  unit  vec- 
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tor  orthogonal  to  i .  Then  the  initial  point 
p3  =(x3,y3,z3)T  of  the  third  segment  is  given  by 


P3  =~d offs*' +( 0,0,  z,  +doff$et  tanyfiT3)T .  (8) 


Impact  Point 


Release  Point 


Wind  direction 


Final  glideslope 


Figure  7  Graphical  representation  of  the  parafoil  trajectory 


Now,  let  pRT  denote  the  position  on  the  RT,  s  de¬ 
note  the  path  length  traveled  along  the  RT  and  VRT  - 
parafoil's  velocity  along  this  path.  Then  the  expression 
for  the  third  segment  of  the  desired  trajectory  is 

pRT(s)  =  p3+T3s,  S=VRT ,  (9) 

where 

T3  =(— cos^f  cosyRV-siny/f  cosyRV-sinyRT3)r  (10) 

is  the  unit  vector  that  points  along  the  line  connecting 
p3  to  pf .  Note,  that  the  length  of  the  third  segment 

s3  =||p3-pf  |  ■  Similarly,  let  s1  denote  the  length  of  the 
first  segment  and  s2  -  of  the  second.  Then,  along  the 
third  segment  ( sl  +  s2  <  s  <  sl  +  s2  +  s3 ) 

pRT(s)  =  p3+T3(s-s1  +  s2) .  (11) 

Next,  consider  Fig. 8.  Let  r  represent  the  radius  of 
the  spiral  descent  segment.  Then  the  center  of  the  spiral 

Pi  =(*i,yi,Zi)T  has  the  following  expression 

p1  =  p3  +  ri*,  (12) 

which  is  used  next  to  determine  the  expression  for  pRT 

along  segment  1.  Let  d  =  iy(x0-x1)2  +  (y0-y1)2 ,  then 


=  sin  1  — ,  Yi  -  tan  l— — — 

d  Xj-x,, 


and  y/2  -  yz-y  +  Ay/  . 


Analogous  to  (10)  define 

Tx  =  ^cos  \// 2  cos^^ ,  si  n  ^  cos^^ , — si  n  Yrt,  j  (13) 
Then 


U 2  _r2 

P2=Po+T1^ - '  S1  =  ||Po  -  P2  II  <14) 

cos  rRTl 

and 

Prt  ($)  =  P0  +TjS  ,  0  <  S  <  S3 .  (15) 

The  turn  rate  y/flr  along  the  spiral  descent  segment 
is  given  by 

Yrt  rt  I"1  cosrRT2-  (16) 


Segment  3  P3 


Figure  8  Top  view  of  the  parafoil  RT 


Then,  along  the  segment  2  ( s1  <s<s1  +  s: 

/  rsin(/7RT (s-sJ  +  a) 
P«r(s)  =  Pi+  -rcos(P;T(s-s1)+A) 

(s1-s)s\nyRT, 


(17) 


where  pKT  —  t,  A  —  i//f+7r  (s2  %)< 


SmX/!T2 

Without  loss  of  generality,  the  z-component  of  the 
center  of  the  spiral  descent  p3  was  set  to  equal  z2  in 

the  definition  of  pflr(s)  above. 

Let 

,  _  A 

dPRT 


T2(s)  = 


ds 


ds 


COs(/?R7- (S  S1)  +  A)cOSyRT2 
sin(/7R*T(s-Si)  +  A)cosyRT2 
-sinyflr, 


.(18) 


v  '  J 

Then,  the  choice  of  A  guarantees  that  the  horizon¬ 
tal  projections  of  vectors  defined  by  Eqns.  (18)  and  (10) 
are  aligned,  i.e. 


rhor 


rhor 


=  (-cosy/f,-siny/f  )T .  (19) 
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This  fact  implies  that  at  p3  the  horizontal  projec¬ 
tions  of  the  commanded  velocity  vectors  for  segments  2 
and  3  are  equal  and  are  independent  of  the  choice  of  r. 

Finally,  the  radius  of  the  spiral  descent  r  is  selected 
to  guarantee  that  at  p2  segments  1  and  2  merge 
smoothly.  This  is  done  numerically  by  solving  a  single 
variable  constrained  optimization  problem.  Note  that  at 
P2 

PfiT  |s=Sl  =P2  =Pi  +  (rsinA,-rcosA,0)T .  (20) 

Therefore,  let 

e(r)  =  ||p2  —  p2  +  (r  si  n  A,  — r  cos  A,  0)r  || ,  (21) 

then  the  desired  value  of  r  is  rdes  =  arg  min e(r ) .  Fig. 9 

rmin  -r-rmax 

includes  a  plot  of  e(r)  versus  r,  which  indicates  that 
there  are  multiple  values  of  rdes .  Therefore  the  values  of 


rmin  and  rmx  are  selected  to  provide  a  unique  solution 


that  is  consistent  with  physical  limitations  of  the  para¬ 
foil.  Note,  in  steady  state  turn  the  bank  angle  of  the 
parafoil  is  defined  by  the  following  expression: 


f Kt  = 


-i  ^ rt_¥rt_  _  tan-i 


Vrt  COS/fi7-2 

rg 


(22) 


where  g  is  acceleration  due  to  gravity.  By  construction 
the  horizontal  projections  of  J1  and  T2  at  p2  are 


aligned,  i.e.  Tr I  =Tf 

1  ls=% 


Figure  9  A  representative  plot  of  e(r)  versus  r 


V.  Integrated  Guidance  and  Control  Algorithm 

The  development  of  the  integrated  guidance  and  control 
algorithm  is  based  on  the  work  reported  in  Ref. 34, 
where  authors  have  proposed  a  new  technique  for  track¬ 
ing  so-called  trimming  trajectories  for  UAV's.  As 


shown  in  Ref. 34  trimming  trajectories  consist  of 
straight  lines  and  helices.  Therefore,  this  methodology 
is  particularly  suitable  for  the  problem  at  hand. 

The  key  ideas  of  the  design  methodology  devel¬ 
oped  in  Ref. 34  include  the  following  five  steps: 

1)  reparameterize  trimming  trajectory  using  the  ar- 
clength  s,  thus  eliminating  time  as  an  independent 
variable; 

2)  resolve  the  position  and  velocity  errors  in  so-called 
F renet  frame  (Ref. 34); 

3)  form  error  dynamics  for  the  system  consisting  of 
the  trajectory  and  parafoil  model,  where  the  posi¬ 
tion  and  velocity  error  states  are  resolved  in  Frenet 
frame; 

4)  design  a  linear  tracking  controller  for  the  lineariza¬ 
tion  of  the  system  along  the  trimming  trajectory; 

5)  implement  this  controller  with  the  true  nonlinear 
plant  using  a  nonlinear  transformation  provided  in 
Ref. 34.  This  implementation  guarantees  that  the 
linearization  along  the  trimming  trajectory  of  the 
feedback  system  consisting  of  the  nonlinear  plant 
and  nonlinear  controller  preserves  the  eigenvalues 
and  transfer  functions  of  the  feedback  interconnec¬ 
tion  of  the  linearized  error  dynamics  and  linear 
tracking  controller  (see  Ref.36). 

V.l  Derivation  of  the  errors  in  Frenet  frame 

As  shown  in  Section  IV  the  desired  RT  is  parameter¬ 
ized  using  the  path/arclength  parameters. 

Let  p  =  (x,y,z)T  denote  current  position  of  the 

parafoil.  The  methodology  developed  in  Ref. 34  requires 
that  position  and  velocity  commands  pfir(s„),  VRT(s,) 
used  to  compute  position  and  velocity  errors, 
pfi7(s,)-p  and  VRT(s,)-V  respectively,  correspond 
to  the  point  on  the  RT  that  is  nearest  to  current  position 
of  the  vehicle.  This  is  done  by  first  determining  the 
value  of  the  path/arclength  parameter 

s,  =  arg  mi n  ||pRr  (s)  —  p||2  and  then  using  it  to  compute 

s 

position  and  velocity  commands.  In  Ref.34  the  problem 
of  determining  s.  is  reduced  to  a  constrained  optimiza¬ 
tion  problem.  However,  computing  limitations  of  the 
onboard  processor  have  imposed  a  need  to  develop  ana¬ 
lytical  techniques  for  the  computation  of  s, ,  discussed 
next. 

An  exact  analytical  expression  for  s,  can  be  de¬ 
rived  for  straight-line  segments  1  and  3.  Recall  Eqns. 
(15)  and  (11).  Then 
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(23) 


s,  =  argmin 


Simple  algebra  shows  that  for  straight-line  seg¬ 
ments 

[T11'(p-Po),0<s<s1, 


^|P0  +T1s-p|, 

0  <  s  <  sx  ' 

Jp3+T3(s-s1-s2)-p||, 

Sx  +  s2  <S<SL  +  S2+S3^ 

(27) 


*-^.y-yi|| 

For  any  vector  pcyJ  =  (xcyl,ycyl,zcyl)T  eC  define  a 


function  n :  R3  ->  R: 


n(pcy()=: 


S.  Ht"  *•  (24) 

T3  (p-p3)  +  (s1  +  s2),s1  +  s2  <s<s1  +  s2  +  s3. 

Now,  recall  that  segment  2  represents  a  spiral.  In 
this  case,  only  an  approximation  of  s,  can  be  found 
analytically.  Its  derivation  is  discussed  next. 

Let 

C  =  {(llm,n)\(l-x1)2  +  (m-y1)2  =r2,z1<n<z2}(25) 

define  a  cylinder  centered  at  p2  =  (x1,y1,z1)T .  Then  the 
spiral  descent  trajectory 

'  rsin(yc>*T(s-s1)  +  A)  ^ 
pfiT(s)  =  p1+  -r  cos(/7fi*T  (s-sJ  +  a)  (26) 

— (s— s1)  sin  ^RT 

V  2  J 

( sx  <  s  <  sx  +  s2 )  is  a  subset  of  C .  F  urthermore,  the  pro¬ 
jection  of  the  vehicle's  position  p  onto  C  is 

r  /  x 

P  proj  ~ 


f  v  \ 

An 

VZny 


, 


rtan" 


■'ey/ 


4, 


Ycyl 


+  7Z 


(28) 


v  cyf  y 

The  function  n  maps  the  cylinder  C  onto  a  rec¬ 
tangle  of  width  2n r  and  height  z2-z1.  M  oreover,  the 
trajectory  pRT  (s)  defined  for  segment  2  is  mapped  into 
a  function  n(pRT(s))  shown  in  Fig. 10. 

Let  d,  denote  the  distance  from  p  to  the  nearest 
point  pRT(s,)on  the  trajectory  pRT (s) .  Then  pRr(s.) 
can  be  obtained  by  determining  the  intersection  of  the 
ball  Bri  of  radius  d,  centered  at  p.  The  mapping 


n(Bc 


d. 

nC  ] 


is  an  ellipse  centered  at  n(p  , )  with 


semi-major  and  semi-minor  radii  given  by  rz  =dxV 3 
r2  +  (r  +  dx)2  -4d2 


and  r=r cos 


2r(r  +  dx 


respectively.  In 


the  previous  expressions  dx  =  ( x — x  . ) 2  +  ( y  —  y 


’  pro i  i 


Intersection  of  this  ellipsoid  with  n(p(s))  represents 
the  mapping  n(p(sj)  of  p(sj  .  Analytical  approxima¬ 
tion  of  n(p(sj)  is  obtained  next  and  then  used  to  com¬ 
pute  s, . 


Figure  10  Graphical  representation  of  n(C ) 


Consider  Fig. 11.  Note  for  the  case  illustrated  in 
this  figure,  the  zn  coordinate  of  n(p(s,))  is  bounded 

above  by  zln  and  below  by  zlp .  Using  the  basic  ge¬ 
ometry  shown  in  Fig. 11,  expressions  for  zln  and  zlp 
are: 


Zln  =  Zproj  +  |mod(z0„  -  zproj ,  2xr  tan  rRT2  )| , 

I  I  (29) 

Zlp~Zproj  ^proj  | COS  / , 

where 
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(31) 


z0n  =z1-s,„sinrRT2 


^*n  ~  Prt 


proj 


■mod(A,2^) 


(30) 


Similarly, 

z2n  =zln-2^rtan  /rTi. 

Z2p  ~  Z  proj  ~  \Z2n  —  Zproj|cOS  /  rj2  • 


The  variables  z2p  and  zlp  represent  the  z- 
coordinates  of  the  projections  of  n(ppfPi )  onto  the  two 
nearest  legs  of  n(p(s)) .  These  variables  can  now  be 
used  to  find  the  approximation  of  zn  ,  zapp : 


app 


^  pip  Zproy|<p2p  Zproj  I’  (32) 

otherwise. 


Note  that  logic  (32)  guarantees  that  in  the  case 
when  two  legs  of  p(s)  are  equidistant  from  p,  the  point 
closer  to  the  ground  is  selected.  Furthermore,  since  n 
mapsz-axisof  R3  onto  itself  zapp  can  be  used  to  com¬ 


pute  s. 


s,  =  — — (33) 
sinrfir2 

which  follows  from  the  definition  of  p(s)  for  the  spiral 
descent  segment.  This  value  of  s,  can  now  be  used  to 
compute  the  unit  basis  vectors  T2,  N2,  and  B2  of  the 
F  renet  frame  {F }  for  segment  2  at  s  =  s, 
dpc  dT2 


ds 


ds 


T2xN2 


*2  _ 

dpc 

,  i^2  - 

dT  2 

’  2  I|t2xn2| 

ds 

ds 

(34) 


or  in  the  final  form 


T,(s.)  == 


v 


COS^y?pp  (s^  S2 )  +  COS^p^ 

sin(Ar(s,-s1)  +  A)cos^T2 
-sinrfiT, 


f 


N,(S.)  = 


sin(A*T(s.-slimi)  +  A) 
cos(/?pp  (s,  —  Sj)  +  a) 

0 


(35) 


T2(5,)xN2(sJ 
2  ||T2(sJx  N2(sJ||' 

For  segments  1  and  3  the  basis  vectors  are  con¬ 
stant.  For  example,  for  segment  1  we  obtain 


^cost//2  cos  rRT^ 

sin 

Ti  = 

si n  y/ 2  cos/RT i 

-  N,= 

COS  y/2 

-sinrfiT 

V  /  RT1  J 

l  o  J 

A  nd  for  segment  3, 


^—COSlf/f  COS/RT^  > 

'  sin^f  ' 

T3  = 

—sin y/ f  cos/R1-3 

-  N3  = 

-cos  y/f 

-sin  vRT 

V  /  rt3  J 

l  0  J 

The  basis  vectors  Bi  and  B3  are  computed  exactly 
as  B2.  Now  the  rotation  matrix  LTFpRi  from  LTP  to  Fre- 
net  frame  for  the  /-th  segment  is  given  by 
ltfpr.  =(Tj,Nf,Bi)7  and  the  corresponding  position 
and  velocity  errors  are 
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Pe  =  trPfRi(P«r(sJ-P)=:(0,ye,Ze)T, 

Pe  =  LTP  (PflT  (^»  )  —  V  )  +  i_jp  Rj  ( Pr7"  (Sj-P). 

N  otice  that  by  construction  the  x-component  of  the  po¬ 
sition  error  vector  pe  is  zero,  see  Ref. 34.  These  error 
vectors  were  used  to  design  a  tracking  control  system 
discussed  in  the  next  section. 

V.2  Control  System  Design 

The  trajectory  tracking  controller  design  methodology 
proposed  in  Ref.34  is  now  applied  to  the  design  of  a 
tracking  controller  for  the  trajectory  generated  using  the 
algorithm  developed  in  Section  IV.  The  controller  can 


only  use  GPS  position  and  velocity  for  feedback.  This 
constraint  is  motivated  by  the  requirement  to  keep  the 
cost  of  the  onboard  avionics  low,  therefore,  only  G  PS  is 
available  for  control. 

During  segments  1  and  2  the  control  system  uses 
differential  flaps  to  drive  the  lateral  components  of  the 
position  and  velocity  error  vectors  to  zero.  On  the  other 
hand,  the  last  segment  is  tracked  using  both  symmetric 
and  differential  flaps,  denoted  here  by  Ss  and  SD ,  re¬ 
spectively.  The  structure  of  tracking  controller  for  the 
lateral  channel  shown  is  Fig. 12. 


The  gains  K  y  and  K  y  shown  in  Fig. 12  were  se¬ 
lected  to  provide  stability  and  performance  for  the  lat¬ 
eral  channel  of  the  feedback  system.  The  variables  z"1 
and  T  denote  backward  shift  and  sampling  period.  The 
sampling  period  for  this  problem  was  dictated  by  the 
onboard  GPS  receiver,  whose  update  rate  is  2Hz. 

The  structure  of  the  vertical  channel  controller  is 
shown  in  Fig. 13.  This  controller  uses  symmetric  flaps 
Ss  to  drive  the  vertical  channel  error  to  0  and  was  only 


used  to  track  the  final  segment  of  the  trajectory.  This 
decision  was  motivated  by  the  fact  that  symmetric  flaps 
have  negligible  control  authority  in  spiral  descent.  Fur¬ 
thermore,  deflecting  symmetric  flaps  tends  to  increase 
the  descent  rate  of  the  parafoil  while  reducing  the  flap 
deflection  budget  available  for  differential  flaps  com¬ 
mand  and,  therefore,  symmetric  flaps  were  not  used 
during  segment  1  as  well. 


Figure  13  Vertical  channel  tracking  controller  structure 


VI.  Simulation  Results 

The  trajectory  tracking  control  system  discussed  in  Sec¬ 
tion  V  has  been  tested  in  simulation  using  a  nonlinear 
6DOF  model  of  the  Pegasus  parafoil  that  also  includes 
the  available  actuator  model.35  Fig. 14  shows  a  compari¬ 


son  between  flight  test  data  collected  during  Pegasus 
drop  at  the  Y  PG  and  model's  response  to  the  wind  col¬ 
lected  during  the  same  drop.  Clearly,  the  6DOF  model 
provides  a  reasonable  approximation  of  the  parafoil 
dynamics. 
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Figure  14  Comparison  of  the  flight  test  data  and  simula¬ 
tion  results 


Fig. 15  shows  results  of  an  extensive  simulation 
analysis  of  the  trajectory  tracking  control  system  using 
the  6DOF  model  discussed  above.  Each  of  40  simula¬ 
tion  runs  included  the  wind  profile  obtained  at  YPG. 
The  GPS  position  and  velocity  errors  were  modeled  as 
white  noise  processes  with  zero  mean  and  standard  de¬ 
viations  shown  in  Tab.l.  Finally,  landing  performance 
statistics  shown  in  Tab.2  indicate  that  system's  per¬ 
formance  exceeded  the  required  circular  error  probable 
(CEP)  of  100m. 


Table  1  GPS  errors  (STD) 


x- direction 

y-di  recti  on 

z-di  recti  on 

Position  (m) 

5 

5 

10 

Velocity  (m/s) 

0.3 

0.3 

0.5 

T able  2  T ouchdown  errors 


CEP (m) 

M  ean  value  (m) 

STD  (m) 

70.7 

67.5 

20.5 

Figure  15  Simulation  results  of  the  trajectory  tracking  control  system 


VII.  Conclusions 

The  optimal  control  strategy  for  the  Pegasus  parafoil 
payload  delivery  system  was  synthesized  based  on 
Pontrjagin's  maximum  principle.  This  motivated  the 
structure  of  the  real-time  reference  trajectory  generator. 
Together  with  the  robust  path  following  algorithm  it 
enabled  a  successful  development  and  simulation  test¬ 
ing  of  the  complete  guidance  and  control  algorithm. 
The  authors  plan  to  flight  test  developed  algorithms  in 
the  early  spring  of  2003  and  to  have  complete  flight  test 
results  soon  afterwards. 
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